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We study the thermophysical properties of warm dense hydrogen using quantum molecular dy- 
namics simulations. New results are presented for the pair distribution functions, the equation of 
state, the Hugoniot curve, and the reflectivity. We compare with available experimental data and 
predictions of the chemical picture. Especially, we discuss the nonmetal-to-metal transition which 
occurs at about 40 GPa in the dense fluid. 
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I. INTRODUCTION 

Hydrogen is an essential element for models of stellar 
and planetary interiors^iS,. The isotopes deuterium and 
tritium are considered as target materials (D-T gas) in 
inertial confinement fusion experiments;^. Therefore, nu- 
merous efforts have been made both experimentally and 
theoretically to understand the behavior of hydrogen, 
deuterium, and tritium in a wide range of densities and 
temperatures. In particular, progress in shock-wave ex- 
perimental technique has allowed the systematic probing 
of the megabar pressure range, so that a sound database 
has been assembled within the last decade. Single or mul- 
tiple shock-wave experiments have been performed for 
hydrogen (or deuterium) by using, e.g., high explosives^, 
gas guns^, pulsed poweri^i^, or high-power lasers^. The 
combination of high pressures and temperatures of sev- 
eral eV defines warm dense matter, a strongly correlated 
state relevant for planetary interiors which is character- 
ized by partial ionization where the bound states exhibit 
a highly transient nature. 

Furthermore, the enormous progress in computer ca- 
pacity has allowed the development and application of 
ab initio simulation techniques for warm dense matter 
such as Path Integral Monte Carlo (PIMC)i^ or Quantum 
Molecular Dynamics (QMD) simulationsii which treat 
quantum effects and correlations systematically. These 
techniques give already highly predictive results for a 
variety of problems and systems; see Ref. [13 for QMD 
simulations. 

The equation of state (EOS) and derived quantities 
such as the Hugoniot curve, the sound velocity, or the 
Griineisen parameter are important material properties 
in this context. Furthermore, optical properties as, e.g., 
the reflectivity are closely related to the dielectric func- 
tion which also determines the dc electrical conductivity 
in the static limit. All these quantities are used to char- 
acterize the unique behavior of warm dense hydrogen, es- 
pecially for high pressures at, or exceeding, one megabar. 



where a transition from a nonconducting, molecular fluid 
to a mono-atomic fluid with metallic-like conductivity 
occurs. Describing the disordered fluid in terms of solid 
state parameters, the fundamental band gap between the 
valence and conduction band decreases with the pres- 
sure and, subsequently, the electrical conductivity shows 
an exponential increase as is typical for thermally acti- 
vated transport in semiconductors"'^'^. For pressures above 
1.4 Mbar, conductivities of about 2000 cm~^, as is 
characteristic for simple metallic fluids such as Cs, have 
been observed experimentally around 3000 K— li^, and 
band gap closure has been claimed to be responsible for 
this nonmetal-to-metal transition. 

On the other hand, concepts of plasma physics have 
been applied to warm dense matter states^^. For in- 
stance, the chemical picture gives a rather simple de- 
scription by identifying stable bound states out of ele- 
mentary particles as new composite particles. Hydrogen 
at normal conditions in this context is a molecular fluid. 
Free electrons are generated at high pressure by dissoci- 
ation of molecules, H2 ^ 2 H, and a subsequent ioniza- 
tion of atoms, H ^ e -I- p. This model yields already 
the strong increase of the conductivity with the pres- 
sure (pressure ionization). In addition, bound states con- 
tribute to conduction via hopping processes^^. The con- 
ceptual problem of all chemical models is the clear defini- 
tion of bound states, the derivation of effective potentials 
between all species, and the calculation of cross sections 
for the respective scattering processes in a strongly cor- 
related medium. 

QMD simulations are a powerful tool to describe warm 
dense matte r^^d^i^°i^^i^^ . The combination of classical 
molecular dynamics for the ions and density functional 
theory (DFT) for the electrons allows one to consider cor- 
relation and quantum effects. Alternatively, wave packet 
simulations have been developed in which the electrons 
are represented on a semi-quantal level by wave packets 

^Typ]\/[p- j23.24.25.26.27.28 

In this paper, we apply QMD simulations and calculate 
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a broad spectrum of thermophysical properties of warm 
dense hydrogen. We determine EOS data for a wide re- 
gion of densities and temperatures and compare with 
chemical models. We calculate the principal Hugoniot 
curve for liquid targets. The Kubo-Greenwood formula 
serves as a starting point for the evaluation of the dy- 
namic conductivity cr(w) from which the dielectric func- 
tion s{uj) and the reflectivity can be extracted. In ad- 
dition, the electronic structure calculation within DFT 
yields the charge density distribution in the simulation 
box at every time step, and the molecular dynamics run 
gives valuable structural information via the ion-ion pair 
correlation function. This is important for the identifi- 
cation and characterization of phase transitions such as 
solid-liquid or liquid-plasma as well as for the nonmetal- 
to-metal transition. 



II. QMD SIMULATIONS 

Within QMD simulations we perform molecular dy- 
namics simulations with a quantum mechanical treat- 
ment of the electrons by using density functional theory 
(DFT). This is based upon the theorems of Hohcnberg 
and Kohn^^ and gives the electron density that mini- 
mizes the ground state energy of the system. It has been 
proven that this density is a unique functional of the ef- 
fective potential VcS- 

From this formalism Kohn and Sham^S derived a com- 
putational scheme which solves the problem for a fictions 
system of non-interacting particles that leads to the same 
electron density. This scheme consists basically of solving 
the Kohn-Sham equations 
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Our ab initio quantum molecular dynamics simula- 
tions were performed within Mermin's finite temperature 
density functional theory (FT-DFT)'^-'^, which is imple- 
mented in the plane wave density functional code VASP 
(Vienna Ab Initio Simulation Package We used 
the projector augmented wave potentials^ and did a gen- 
eralized gradient approximation (GGA) using the param- 
eterization of PBE^i. Extensive test calculations, as per- 
formed already by Desjarlais^^, have shown that the EOS 
data are dependent on the plane wave cutoff. A conver- 
gence of better than 1% is secured for Ecut — 1200 eV 
which was used in all actual calculations. The elec- 
tronic structure calculations were performed for a given 
array of ion positions which are subsequently varied by 
the forces obtained within the DFT calculations via the 
Hellmann-Feynman theorem for each molecular dynam- 
ics step. This schema is repeated until the EOS mea- 
sures are converged and a thermodynamic equilibrium is 
reached. 



The simulations were done for 64 atoms in a supercell 
with periodic boundary conditions. The temperature of 
the ions was controled by a Nose thermostalj^ and the 
temperature of the electrons was fixed by Fermi weighting 
the occupation of bands'^'^. The Brillouin zone was sam- 
pled by evaluating the results at Baldereschi's mean value 
point"^^ which showed best agreement with a sampling of 
the Brillouin zone using a higher number of k-points. 
The density of the system was fixed by the size of the 
simulated supercell. To achieve a small statistical error 
due to fluctuations the system was simulated 1000-1500 
steps further after reaching the thermodynamic equilib- 
rium. The EOS data and pair correlation functions were 
then obtained by averaging over all particles and sim- 
ulation steps in equilibrium. Similar calculations were 
performed recently for the thermophysical properties of 
warm dense helium^S in order to verify the nonmetal-to- 
metal transition at high pressures. 

The zero-point vibrational energy of the H2 molecules 
is not included in DFT calculations. In previous calcu- 
lations, the energy ^hvmb per molecule is simply added 
which is very important, especially at low temperatures 
and for the calculation of an exact initial internal energy 
for the reference state of the Hugoniot curve, which is 
0.0855 g/cvc? at 20 K. To account for this quantum effect 
more sensitively for arbitrary temperatures, the fraction 
of molecules has to be derived, e.g., for all states along 
the Hugoniot curve. This can be done via the coordina- 
tion number 



K{r) = 



N -I 
V 



Attt'"^ g{r')dr' , 



(2) 



which is a weighted integral over the pair correlation 
function g{r) of the ions. N denotes the number of ions 
and V the volume of the supercell in the simulation. The 
doubled value of K at the maximum of the molecular 
peak in g{r)^ which is found around r = 0.748 A, is 
then equal to the fraction of ions bound to a molecule 
and twice the amount of molecules in the supercell. An 
example is shown in Fig. [1] where the increasing dissoci- 
ation with higher density can be seen. In Fig. [2] we show 
the thermal dissociation; the molecular peak dissappears 
with increasing temperature at constant density. Note 
that the peak is thermally broadened. 

The dissociation degree is calculated for a number of 
isotherms and then approximated by a Fermi function 
which has two adjustable parameters. These parame- 
ters can be represented by temperature-dependent func- 
tions so that the dissociation degree and, subsequently, 
the contribution of molecules to the zero-point inter- 
nal energy are determined for arbitrary temperatures. 
The results show that molecules can be neglected above 
10,000 K. 

We compare the resulting dissociation degree with that 
derived by Vorberger et al^^ in Fig. [3l They counted all 
pairs of atoms in a range of 1.8 a^ as atoms. In a second 
step they reduced the number of molecules by counting 
only those pairs that are stable for longer than ten vibra- 
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FIG. 1: Proton-proton pair distribution function and cor- 
responding coordination numbers according to Eq. ([2} for 
1000 K and three densities. 
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FIG. 2: Proton-proton pair distribution function and cor- 
responding coordination numbers according to Eq. ([2} for 
0.5 g/cm'^ and three temperatures. 



tional periods. In all three cases the amount of molecules 
is lower for higher densities und the molecules disap- 
pear at higher temperatures due to thermal dissociation. 
This picture shows that the dissociation degree depends 
strongly on the definition of the term molecule in the 
warm dense matter region. Our alternative method gives 
a smoother behavior of the dissociation degree which 
starts at lower temperatures and is in between the two 
cases described by Vorberger ei al^ at higher tempera- 
tures. 



III. RESULTS FOR THE EOS AND HUGONIOT 
CURVES 

We show the thermal EOS of warm dense hydrogen in 
Fig. m The isotherms of the pressure show a systematic 
behavior in terms of the density and temperature. We 
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FIG. 3: Ratio of hydrogen molecules with respect to the to- 
tal number of protons for three densities: Our coordination 
number method (solid) is compared with the pair-counting 
method of Vorberger et al.'^'^ (dotted). Their result count- 
ing only pairs with a lifetime longer than ten H2 vibrational 
periods is also given (dashed line). 



find no indication of a plasma phase transition (PPT) 
which would result in an instability of the EOS isotherms 
that would need to be treated by a Maxwell construc- 
tion^. The absence of a PPT is in contrast to results of 
chemical models which use, e.g., Fluid Variational The- 
ory (FVT)^i^i^ or liquid state perturbation theory^. 
Chemical models are based on a free energy minimization 
schema for a mixture of hydrogen atoms, molecules, and 
a plasma in chemical equilibrium. Correlations are taken 
into account based on effective two-particle potentials. 
The description of the free charged particles (plasma) is 
done beyond the Debye-Hiickel approximation by using 
efficient Pade formulas which are valid for a wide region 
of densities and temperatures. 

In Fig.|3]our QMD results are compared with the chem- 
ical models FVTii and SCvH-i^. The EOS derived by 
Saumon et al^ shows also a PPT (SCvH-ppt data set). 
The modified SCvH-i data set shown here avoids the PPT 
by using an interpolation through the instability region. 
Therefore, both data sets can be used to study the influ- 
ence of a PPT on interior models of giant planets such as 
Jupiter. Consistent chemical models yield the correct 
low-temperature and low-density limit and agree with 
our QMD results there. A good agreement is also found 
in the high-density limit where a nearly temperature in- 
dependent behavior characteristic of a degenerate plasma 
is found. At medium densities the pressure isotherms of 
FVT and SCvH-i he well below the QMD data; the de- 
viations amount up to 25%. 

We have encountered a region with {dP/dT)v < 0, 
which was previously reported by Vorberger et al^. It 
can be related to the rapid dissociation transition at low 
temperatures. 

Following Lenosky et al^ and Beule et al^ we fit 
smooth functions for the pressure P and the internal en- 
ergy U as an expansion in terms of density p and tem- 
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TABLE L Coefficients a^fc in the expansion for the pressure 
P'"' according to Eqs. (g)) and (O. 
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TABLE II: Coefficients bjk in the expansion for the specific 
internal energy u according to Eqs. © and ([7|). 
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perature T to the given results of the QMD simulations. 
The pressure is split into an ideal and an interaction con- 
tribution: 

p ^ p^d ^ pM ^ ^ p^nt / ^ rpy r^. 

rriH 

The QMD data for the pressure P given in kbar can be 
interpolated by the following expansion for the interac- 
tion contribution: 

P-*(p, T) - (Ai(r) + ^2(T)p)^«(^) , (4) 



Ai{T) = Ojo cxp ^- ^ "^^ "'^ ^ j + a,j3 + auT. (5) 

The coefficients a^fc are summarized in Tab. [D 

The QMD data for the specific internal energy u = 
U/m given in kJ/g can be interpolated by a similar ex- 
pansion: 



(6) 



Bj{T) = bjo exp 



b,2 . 



+ 6,3 + h^T. (7) 



The expansion coefficients are given in Tab. [TTl 

The expansions and ((S]) reproduce the ah initio 
QMD data within 5% accuracy in a density range from 
0.5 g/cm^ to 5 gjcvc? between 500 K and 20000 K and can 
easily be applied in planetary models or hydrodynamic 
simulations for warm dense matter. The expansions fulfill 
thermodynamic consistency expressed by the relation 



P -T 



dP 
df 



dV 



(8) 




P [g/cm ] 



FIG. 4: Thermal EOS for warm dense hydrogen (pressure 
isotherms): QMD data are compared with the chemical mod- 
els FVT+ii and SCvH-i*^ 



within 15% accuracy which is mainly due to the devia- 
tions from the QMD data itself. 

A crucial measure for theoretical EOS data is the prin- 
cipal Hugoniot curve which is plotted in Fig. [5l It de- 
scribes all possible final states (p, P, u) of shock wave ex- 
periments according to the Hugoniot equation 



2 Po Po 



(9) 



starting at the same initial conditions {pq,Pq,uq). For 
the hydrogen principal Hugoniot curve, the initial density 
is Po — 0.0855 g/cm'^ and the initial internal energy uq — 
—314 kJ/g at a temperature of 20 K. The initial pressure 
Po can be neglected because of the high pressure of the 
final state. 

Shock wave experiments have been performed for deu- 
terium using gas guns^, magnetically launched flyer 
plates at Sandia's Z machine^ or high explosives (HE)i^. 
These experiments indicate a maximum compression of 
4.25 at about 50 GPa. 

Another series of laser-driven experiments^ shows sys- 
tematic deviations from the experiments quoted above. 
Especially, a maximum compression of 6 has been re- 
ported at about 1 Mbar. According to the unani- 
mous evaluation of the shock-wave experimental data for 
molecular liquids'' , we compare our QMD data in Fig. O 
only with the data sets mentioned above. 

The systematic increase of the cutoff energy Ecut in 
QMD simulations from 500 eV^ to 1200 eV^ has lead 
to fully converged results in agreement with the exper- 
imental points. The consideration of the zero-point vi- 
brations of the H2 molecules along the entire Hugoniot 
curve yields a very good agreement of QMD data with 
the gas gun experiments^ expecially for low pressures. 
The calculated Hugoniot curve has a maximum compres- 
sion of 4.5 which is slightly higher than the HE and Z 
experiments indicate (about 4.25). This is an agreement 
of about 5% accuracy which can be translated into an ac- 
curacy of about 1% in the measured shock and particle 
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velocity, which is in the range of the systematic errors in 
the experiments. The compression decreases with higher 
pressures and temperatures and reaches the correct high- 
temperature hmit as given by the PIMC simulations^^. 
The QMD curve lies slightly below the experimental data 
for compression rates between 3 and 4 which could be due 
to the known band gap problem of DFT in GGA. The 
FVT curve**^ is shown as a representative of chemical 
models which, in general, show a higher compressibility 
well beyond 4.5. 

Also shown is the linear mixing result of Ross^. This 
curve shows a sixfold compression and is not in agreement 
with the shown experiments. The curve of Kerlej*^ has 
a maximum compression of 4.25, like the experiments 
indicate, but the pressure is there slightly higher than 
the results of the QMD simulations. 



200 




FIG. 5: Principal Hugoniot curve for hydrogen. The re- 
sults of this work (solid line) are compared with previous 
QMD results of Lenosky et al^ (dashed) and Desjarlaia^ 
(stars), PIMC simulations^ (dotted), the linear mixing model 
of Ross^ (dot-dash-dashed), the model of Kerley^ (dot-dot- 
dashed) and the chemical model FVT'^^ (dot-dashed). Ex- 
periments: Gas gun** (diamonds), Sandia Z machine* (grey 
squares; grey line: running average through the Us-Up data), 
high explosives^ (black circles). 



IV. DYNAMIC CONDUCTIVITY, 
REFLECTIVITY AND DC CONDUCTIVITY 

The dynamic conductivity (jij-ij) is derived from the 
Kubo-Greenwood formula i^^i^^ 



cr(w) 



N N 3 

E ^c^) E E E [^(^-k) - ne,,k)] 

k i — 1 i—1 cx—1 



x|(*j,k|Va|*j^k)|^(5(ej,k - e»,k - f>w), 



(10) 



where e is the electron charge and m its mass. The sum- 
mations over i and j run over N descrete bands consid- 
ered in the electronic structure calculation for the cubic 
supercell volume fl. The three spatial directions are av- 
eraged by the a sum. F{ei^k) describes the occupation 



of the ith band corresponding to the energy e^^k and the 
wavefunction 'I'i^k at k. The i5-function has to be broad- 
ened because a discrete energy spectrum results from 
the finite simulation volume^'^. Integration over the Bril- 
louin zone is performed by sampling special k points'"'^, 
where Vl^(k) is the respective weighting factor. We used 
Baldereschi's mean value poin1^ to reach a convergence 
of better than 10% accuracy. 

Optical properties can be derived from the frequency- 
dependent conductivity Eq. (llOp . The standard method 
is to obtain the imaginary part via the Kramers-Kronig 
relation 



TT J [V^-U^) 



(11) 



P is the prinicipal value of the integral. The dielectric 
function can be calculated directly with the conductivity: 



ei(w) = 1 3_cr2(tj), 



(12) 
(13) 



The square of the index of refraction contains the real 
part n and the imaginary part k is equal to the dielectric 
function which leads to the following relations: 



1 



= -v/|e(^)| + |ei(a;)|, 



k{Lo) = -VkHI-kiMI. 



(14) 
(15) 



The index of refraction is then used to calculate optical 
propersties such as the reflectivity r: 



r(w) 



[I - n{Lo)f + kjujf 
[l + n{uj)Y^k{uf 



(16) 



We compare our ah initio results with reflectivities 
measured along the Hugoniot curve5§- in Fig. [SI the agree- 
ment is excellent. The change of the hydrogen reflectivity 
with the pressure can be interpreted as a gradual transi- 
tion from a molecular insulating fluid through an atomic 
fluid above 20 GPa where the atoms have strongly fluc- 
tuating bonds with next neighbors^° to a dense, almost 
fully ionized plasma with a reflectivity of about 50-60 % 
at high pressures above 40 GPa. The chemical model- 
shows also this qualitative behavior but the abrupt in- 
crease of the reflectivity occurs at a higher density. This 
shows the difficulties of the chemical models in finding 
the correct shifts of the dissociation and ionization en- 
ergies as function of density and temperature and, thus, 
the location of the nonmetal-to-metal transition. How- 
ever, the limits of a molecular fluid at low pressures and 
of a fully ionized plasma at high pressures are incorpo- 
rated in a reasonable way. 
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FIG. 6: Reflectivity for a wavelength of 808 nm along the 
Hugoniot curve of hydrogen: QMD results are compared with 
experimental data of Celliers et al^ and predictions of the 
chemical model FVT''^ using the COMPTRA code^. 



V. CONCLUSION 

We have performed ah initio QMD simulations to 
study thermophysical properties of hydrogen under ex- 
treme conditions. As a result we obtained highly con- 
verged EOS data which are relevant for modeling giant 
planets and for the understanding of the fundamental 
behavior of hydrogen at high pressure. The deviations 
between our QMD data and chemical models amount up 
to 25%. We have constructed smooth fit functions for the 
QMD data for the pressure and the internal energy which 



can be used easily in, e.g., hydrodynamic simulations for 
warm dense hydrogen and in astrophysical applications. 

The results show a smooth transition from a molecular 
liquid to an atomic fluid of metal-like state. There were 
no signs of a PPT which is predicted by other models. 

With these EOS results we have calculated the princi- 
pal Hugoniot curve which is in agreement with dynamic 
experiments and has the correct high-temperature limit 
as given by PIMC simulations. 

We obtained optical properties using the Kubo- 
Greenwood formula. The reflectivity along the Hugoniot 
curve is in excellent agreement with experiments. The 
results show the occurrence of a nonmetal-to-metal tran- 
sition at about 40 GPa. 
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